Load packages

library(ggplot2)
library(cowplot)
# library(lme4)
suppressWarnings(suppressMessages(library("lme4")))
library(pbkrtest)
library(tidyr)
library(dplyr)
library(optimx)
library(broom)
library(vegan)

Data

herbdat <- read.csv("data/Herb_20170303.csv")

summary(herbdat)
##       site            code          quadrat     
##  Min.   :1.000   Min.   :101.0   Min.   :101.0  
##  1st Qu.:1.000   1st Qu.:102.0   1st Qu.:203.0  
##  Median :2.000   Median :104.0   Median :305.0  
##  Mean   :2.003   Mean   :103.9   Mean   :300.9  
##  3rd Qu.:3.000   3rd Qu.:106.0   3rd Qu.:408.0  
##  Max.   :3.000   Max.   :107.0   Max.   :510.0  
##                                                 
##                                 species        no.ind.       
##  Cardamine leucantha                :1084   Min.   :  1.000  
##  Meehania urticifolia               : 817   1st Qu.:  2.000  
##  Brachybotrys paridiformis          : 543   Median :  4.000  
##  Osmorhiza aristata                 : 538   Mean   :  9.718  
##  Oxalis acetosella subsp. Griffithii: 440   3rd Qu.: 10.000  
##  Pseudostellaria sylvatica          : 427   Max.   :276.000  
##  (Other)                            :6547                    
##      height           coverage       census     treatment   exclosure    
##  Min.   :  1.000   Min.   :0.0000   15fa:1295   C:2228    Min.   :0.000  
##  1st Qu.:  4.000   1st Qu.:0.0500   16fa:2083   F:2123    1st Qu.:0.000  
##  Median :  6.000   Median :0.1300   16sp:4183   I:2182    Median :1.000  
##  Mean   :  8.881   Mean   :0.2007   16su:2835   S:1634    Mean   :0.504  
##  3rd Qu.: 11.000   3rd Qu.:0.3200               W:2229    3rd Qu.:1.000  
##  Max.   :126.000   Max.   :0.9100                         Max.   :1.000  
##  NA's   :4         NA's   :132
unique(herbdat$species)
##  [1] Filipendula palmata                
##  [2] Cardamine leucantha                
##  [3] Meehania urticifolia               
##  [4] Osmorhiza aristata                 
##  [5] Athyrium brevifrons                
##  [6] Oxalis acetosella subsp. Griffithii
##  [7] Circaea cordata                    
##  [8] Lamium barbatum                    
##  [9] Impatiens noli-tangere             
## [10] Carex pilosa                       
## [11] Geum aleppicum                     
## [12] Arisaema amurense                  
## [13] Viola verecunda                    
## [14] Plagiorhegma dubia                 
## [15] Galium triflorum                   
## [16] Diarrhena manshurica               
## [17] Brachybotrys paridiformis          
## [18] Equisetum hyemale                  
## [19] Circaea lutetiana                  
## [20] Osmunda cinnamomea                 
## [21] Parasenecio hastatus               
## [22] Paeonia obovata                    
## [23] Phryma leptostachya subsp. asiatica
## [24] Saussurea manshurica               
## [25] Paris verticillata                 
## [26] Adiantum pedatum                   
## [27] Dryopteris crassirhizoma           
## [28] Convallaria majalis                
## [29] unknown spp1                       
## [30] Carex callitrichos                 
## [31] Carex siderosticta                 
## [32] Dioscorea nipponica                
## [33] Angelica amurensis                 
## [34] Carex remotiuscula                 
## [35] unknown spp2                       
## [36] Adenophora remotiflora             
## [37] Artemisia stolonifera              
## [38] Pseudostellaria sylvatica          
## [39] Maianthemum bifolium               
## [40] Polemonium coeruleum               
## [41] Viola acuminata                    
## [42] Matteuccia struthiopteris          
## [43] Anemone baicalensis                
## [44] Chrysosplenium sinicum             
## [45] Aruncus sylvester                  
## [46] Doellingeria scaber                
## [47] Sanicula rubriflora                
## [48] Rubia cordifolia                   
## [49] Smilacina japonica                 
## [50] Actaea asiatica                    
## [51] Allium monanthum                   
## [52] Anemone spp.                       
## [53] Adoxa moschatellina                
## [54] Corydalis repens                   
## [55] Anemone amurensis                  
## [56] Enemion raddeanum                  
## [57] Corydalis turtschaninovii          
## [58] Lilium distichum                   
## [59] Anemone raddeana                   
## [60] Thalictrum spp.                    
## [61] Trillium kamtschaticum             
## [62] Anemone reflexa                    
## [63] Ranunculus japonicus               
## [64] Astilbe chinensis                  
## [65] Veratrum nigrum                    
## [66] Adonis amurensis                   
## [67] unknown spp4                       
## [68] Gagea triflora                     
## [69] Gagea pauciflora                   
## [70] Eranthis stellata                  
## [71] Aconitum kusnezoffii               
## [72] Aconitum sczukinii                 
## [73] Hylomecon japonica                 
## [74] Corydalis yanhusuo                 
## [75] Urtica laetevirens                 
## [76] Adenocaulon himalaicum             
## [77] Rabdosia excisa                    
## [78] unknown spp5                       
## [79] Chrysosplenium lectus-cochleae     
## [80] Pilea pumila                       
## [81] Agrimonia pilosa                   
## [82] unknown spp3                       
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
hist(table(herbdat$species))

## manipulate data
herbdat$site <- factor(herbdat$site)
herbdat$quadrat <- factor(herbdat$quadrat)
herbdat$census <- factor(herbdat$census)
herbdat$exclosure <- factor(herbdat$exclosure)

## Removing this to prevent problems later.
herbdat$species <- factor(trimws(as.character(herbdat$species)))

herbdat$quad.unique <- as.factor(paste(herbdat$site,
                                          herbdat$exclosure,
                                          herbdat$quadrat, sep='_'))
herbdat$plot <- as.factor(paste(herbdat$site, herbdat$exclosure, 
                                   sep='_'))
## calculate diversity
sp.diversity_herb <- summarise(
  group_by(herbdat, plot, quadrat, quad.unique, treatment, exclosure, census), 
  shannon=diversity(no.ind.), simpson=diversity(no.ind., index='invsimpson'))
sp.diversity_herb
## Source: local data frame [1,190 x 8]
## Groups: plot, quadrat, quad.unique, treatment, exclosure [?]
## 
##      plot quadrat quad.unique treatment exclosure census   shannon
##    <fctr>  <fctr>      <fctr>    <fctr>    <fctr> <fctr>     <dbl>
## 1     1_0     101     1_0_101         W         0   15fa 1.1988493
## 2     1_0     101     1_0_101         W         0   16fa 1.6512596
## 3     1_0     101     1_0_101         W         0   16sp 1.3934444
## 4     1_0     101     1_0_101         W         0   16su 1.4727372
## 5     1_0     102     1_0_102         C         0   15fa 0.5623351
## 6     1_0     102     1_0_102         C         0   16fa 1.0750931
## 7     1_0     102     1_0_102         C         0   16sp 1.4268571
## 8     1_0     102     1_0_102         C         0   16su 1.3457086
## 9     1_0     103     1_0_103         I         0   15fa 1.0296530
## 10    1_0     103     1_0_103         I         0   16fa 1.2530462
## # ... with 1,180 more rows, and 1 more variables: simpson <dbl>
sp.diversity_herb$census <- factor(sp.diversity_herb$census,
                                   levels=c('15fa','16sp','16su','16fa'))


## data set of pesticide and snow removal treatment seperately
## pesticide
herbpest <- subset(herbdat, treatment !='S')
herbpest$treatment <- factor(herbpest$treatment, levels=c('W', 'F', 'I', 'C'))

sp.diversity_herbpest <- subset(sp.diversity_herb, treatment !='S')
sp.diversity_herbpest$treatment <- factor(sp.diversity_herbpest$treatment, levels=c('W', 'F', 'I', 'C'))
summary(sp.diversity_herbpest)
##   plot        quadrat     quad.unique  treatment exclosure  census   
##  1_0:152   209    : 25   1_1_209:  5   W:237     0:472     15fa:240  
##  1_1:161   504    : 25   2_1_504:  5   F:237     1:482     16sp:241  
##  2_0:160   101    : 24   1_0_101:  4   I:240               16su:233  
##  2_1:161   102    : 24   1_0_102:  4   C:240               16fa:240  
##  3_0:160   103    : 24   1_0_103:  4                                 
##  3_1:160   104    : 24   1_0_104:  4                                 
##            (Other):808   (Other):928                                 
##     shannon         simpson      
##  Min.   :0.000   Min.   : 1.000  
##  1st Qu.:1.278   1st Qu.: 2.794  
##  Median :1.660   Median : 3.984  
##  Mean   :1.603   Mean   : 4.353  
##  3rd Qu.:1.967   3rd Qu.: 5.490  
##  Max.   :2.973   Max.   :15.979  
## 
## snow
herbsnow <- subset(herbdat, treatment =='C' | treatment =='S')
sp.diversity_herbsnow <- subset(sp.diversity_herb, treatment =='C' | treatment =='S')

1. Analysis of diversity for pesticide

ggplot(sp.diversity_herbpest, aes(x=treatment, y=shannon, colour = exclosure)) +
  geom_boxplot()

summarise(group_by(sp.diversity_herbpest, treatment, exclosure, census), meanH=mean(shannon), 
          seH=sd(shannon)/sqrt(length(shannon))) %>%
  ggplot(aes(x=treatment, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census, 
             shape=exclosure)) + geom_pointrange() + coord_trans(y="exp") +
  labs(x="pesticide", y="Effective number of species")

mod.pest.herb.div1 <- lmer(shannon ~ census*treatment*exclosure + (1|plot/quad.unique), data=sp.diversity_herbpest)
mod.pest.herb.div2 <- lmer(shannon ~ census*(treatment + exclosure) + (1|plot/quad.unique), data=sp.diversity_herbpest)
mod.pest.herb.div3 <- lmer(shannon ~ census + treatment*exclosure + (1|plot/quad.unique), data=sp.diversity_herbpest)
mod.pest.herb.div4 <- lmer(shannon ~ census*treatment + exclosure + (1|plot/quad.unique), data=sp.diversity_herbpest)
mod.pest.herb.div5 <- lmer(shannon ~ census + treatment + exclosure + (1|plot/quad.unique), data=sp.diversity_herbpest)
anova(mod.pest.herb.div1, mod.pest.herb.div2, mod.pest.herb.div3,
      mod.pest.herb.div4, mod.pest.herb.div5)
## refitting model(s) with ML (instead of REML)
## Data: sp.diversity_herbpest
## Models:
## mod.pest.herb.div5: shannon ~ census + treatment + exclosure + (1 | plot/quad.unique)
## mod.pest.herb.div3: shannon ~ census + treatment * exclosure + (1 | plot/quad.unique)
## mod.pest.herb.div4: shannon ~ census * treatment + exclosure + (1 | plot/quad.unique)
## mod.pest.herb.div2: shannon ~ census * (treatment + exclosure) + (1 | plot/quad.unique)
## mod.pest.herb.div1: shannon ~ census * treatment * exclosure + (1 | plot/quad.unique)
##                    Df    AIC     BIC  logLik deviance   Chisq Chi Df
## mod.pest.herb.div5 11 854.15  907.61 -416.07   832.15               
## mod.pest.herb.div3 14 859.31  927.36 -415.65   831.31  0.8379      3
## mod.pest.herb.div4 20 861.88  959.09 -410.94   821.88  9.4266      6
## mod.pest.herb.div2 23 832.69  944.49 -393.35   786.69 35.1884      3
## mod.pest.herb.div1 35 846.64 1016.77 -388.32   776.64 10.0479     12
##                    Pr(>Chisq)    
## mod.pest.herb.div5               
## mod.pest.herb.div3     0.8404    
## mod.pest.herb.div4     0.1510    
## mod.pest.herb.div2  1.112e-07 ***
## mod.pest.herb.div1     0.6118    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod.pest.herb.div2, type=c('p', 'smooth'))

dotplot(ranef(mod.pest.herb.div2, condVar=TRUE))
## $`quad.unique:plot`

## 
## $plot

library(sjPlot)
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
sjp.lmer(mod.pest.herb.div2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.pest.herb.div2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.pest.herb.div2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## shannon ~ census * (treatment + exclosure) + (1 | plot/quad.unique)
##    Data: sp.diversity_herbpest
## 
## REML criterion at convergence: 867.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7833 -0.5039  0.0846  0.5760  2.6205 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.06172  0.2484  
##  plot             (Intercept) 0.02127  0.1458  
##  Residual                     0.09811  0.3132  
## Number of obs: 954, groups:  quad.unique:plot, 240; plot, 6
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)            1.2862267  0.1020736  12.601
## census16sp             0.5759204  0.0639256   9.009
## census16su             0.6278072  0.0653146   9.612
## census16fa             0.3230473  0.0639365   5.053
## treatmentF            -0.0197691  0.0728158  -0.271
## treatmentI            -0.0472976  0.0729310  -0.649
## treatmentC            -0.1203521  0.0728752  -1.651
## exclosure1            -0.1504535  0.1297775  -1.159
## census16sp:treatmentF  0.0490280  0.0808740   0.606
## census16su:treatmentF -0.1126312  0.0821542  -1.371
## census16fa:treatmentF -0.0464141  0.0808740  -0.574
## census16sp:treatmentI  0.0233930  0.0807355   0.290
## census16su:treatmentI -0.0231689  0.0817315  -0.283
## census16fa:treatmentI -0.0603406  0.0808740  -0.746
## census16sp:treatmentC  0.1706872  0.0808740   2.111
## census16su:treatmentC  0.0002213  0.0816560   0.003
## census16fa:treatmentC  0.0573980  0.0808740   0.710
## census16sp:exclosure1  0.3074598  0.0571376   5.381
## census16su:exclosure1  0.1070975  0.0577662   1.854
## census16fa:exclosure1  0.0291159  0.0571865   0.509
## 
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

2. Analysis of diversity for snow removal treatment

ggplot(sp.diversity_herbsnow, aes(x=treatment, y=shannon, colour = exclosure)) +
  geom_boxplot()

summarise(group_by(sp.diversity_herbsnow, treatment, exclosure, census), meanH=mean(shannon), 
          seH=sd(shannon)/sqrt(length(shannon))) %>%
  ggplot(aes(x=treatment, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census, 
             shape=exclosure)) + geom_pointrange() + coord_trans(y="exp") +
  labs(x="pesticide", y="Effective number of species")

mod.snow.herb.div1 <- lmer(shannon ~ census*treatment*exclosure + (1|plot/quad.unique), data=sp.diversity_herbsnow)
mod.snow.herb.div2 <- lmer(shannon ~ census*(treatment + exclosure) + (1|plot/quad.unique), data=sp.diversity_herbsnow)
mod.snow.herb.div3 <- lmer(shannon ~ census + treatment*exclosure + (1|plot/quad.unique), data=sp.diversity_herbsnow)
mod.snow.herb.div4 <- lmer(shannon ~ census*treatment + exclosure + (1|plot/quad.unique), data=sp.diversity_herbsnow)
mod.snow.herb.div5 <- lmer(shannon ~ census + treatment + exclosure + (1|plot/quad.unique), data=sp.diversity_herbsnow)
anova(mod.snow.herb.div1, mod.snow.herb.div2, mod.snow.herb.div3,
      mod.snow.herb.div4, mod.snow.herb.div5)
## refitting model(s) with ML (instead of REML)
## Data: sp.diversity_herbsnow
## Models:
## mod.snow.herb.div5: shannon ~ census + treatment + exclosure + (1 | plot/quad.unique)
## mod.snow.herb.div3: shannon ~ census + treatment * exclosure + (1 | plot/quad.unique)
## mod.snow.herb.div4: shannon ~ census * treatment + exclosure + (1 | plot/quad.unique)
## mod.snow.herb.div2: shannon ~ census * (treatment + exclosure) + (1 | plot/quad.unique)
## mod.snow.herb.div1: shannon ~ census * treatment * exclosure + (1 | plot/quad.unique)
##                    Df    AIC    BIC  logLik deviance   Chisq Chi Df
## mod.snow.herb.div5  9 454.31 491.80 -218.15   436.31               
## mod.snow.herb.div3 10 455.68 497.34 -217.84   435.68  0.6253      1
## mod.snow.herb.div4 12 451.36 501.34 -213.68   427.36  8.3273      2
## mod.snow.herb.div2 15 442.12 504.60 -206.06   412.12 15.2329      3
## mod.snow.herb.div1 19 447.76 526.90 -204.88   409.76  2.3668      4
##                    Pr(>Chisq)   
## mod.snow.herb.div5              
## mod.snow.herb.div3   0.429083   
## mod.snow.herb.div4   0.015551 * 
## mod.snow.herb.div2   0.001628 **
## mod.snow.herb.div1   0.668630   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod.snow.herb.div2, type=c('p', 'smooth'))

dotplot(ranef(mod.snow.herb.div2, condVar=TRUE))
## $`quad.unique:plot`

## 
## $plot

sjp.lmer(mod.snow.herb.div2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.snow.herb.div2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.snow.herb.div2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## shannon ~ census * (treatment + exclosure) + (1 | plot/quad.unique)
##    Data: sp.diversity_herbsnow
## 
## REML criterion at convergence: 456.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.05453 -0.54855  0.03804  0.62872  2.02116 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.09042  0.3007  
##  plot             (Intercept) 0.03795  0.1948  
##  Residual                     0.09342  0.3057  
## Number of obs: 476, groups:  quad.unique:plot, 122; plot, 6
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            1.18462    0.13124   9.027
## census16sp             0.78720    0.06853  11.487
## census16su             0.66114    0.06949   9.514
## census16fa             0.42863    0.06853   6.255
## treatmentS            -0.29784    0.07830  -3.804
## exclosure1            -0.17043    0.17750  -0.960
## census16sp:treatmentS -0.08247    0.07957  -1.036
## census16su:treatmentS  0.13544    0.08011   1.691
## census16fa:treatmentS  0.09719    0.07957   1.221
## census16sp:exclosure1  0.22628    0.07955   2.844
## census16su:exclosure1  0.05764    0.08012   0.719
## census16fa:exclosure1 -0.06725    0.07955  -0.845

3. Abundance analysis

First need to calculate abundance. In this analysis, we examine for overall, pesticide and snow removal treatment seperately, and then move to individual species.

## Overall
## data set
herb.abundat <- aggregate(no.ind. ~ treatment + exclosure + site + quadrat + plot + quad.unique + census, data=herbdat, FUN=sum)
herb.abundat$no.ind. <- as.numeric(herb.abundat$no.ind.)

mod.herb.abun1 <- glmer(no.ind. ~ site + treatment*exclosure + census +(1|quad.unique),
                       data=herb.abundat, family=poisson, 
                       control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.herb.abun2 <- glmer(no.ind. ~ site + treatment*exclosure + census +(1|plot/quad.unique),
                        data=herb.abundat, family=poisson,
                        control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.herb.abun3 <- glmer(no.ind. ~ site + census*(treatment + exclosure) +(1|plot/quad.unique),
                        data=herb.abundat, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.herb.abun4 <- glmer(no.ind. ~ site + census + treatment*exclosure +(1|plot/quad.unique),
                        data=herb.abundat, family=poisson,
                        control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.herb.abun5 <- glmer(no.ind. ~ site + census + treatment + exclosure +(1|plot/quad.unique),
                        data=herb.abundat, family=poisson,
                        control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.herb.abun1, mod.herb.abun2, mod.herb.abun3,
      mod.herb.abun4, mod.herb.abun5)
## Data: herb.abundat
## Models:
## mod.herb.abun5: no.ind. ~ site + census + treatment + exclosure + (1 | plot/quad.unique)
## mod.herb.abun1: no.ind. ~ site + treatment * exclosure + census + (1 | quad.unique)
## mod.herb.abun2: no.ind. ~ site + treatment * exclosure + census + (1 | plot/quad.unique)
## mod.herb.abun4: no.ind. ~ site + census + treatment * exclosure + (1 | plot/quad.unique)
## mod.herb.abun3: no.ind. ~ site + census * (treatment + exclosure) + (1 | plot/quad.unique)
##                Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.herb.abun5 13 15060 15126 -7516.7    15034                         
## mod.herb.abun1 16 15074 15156 -7521.2    15042  0.000      3          1
## mod.herb.abun2 17 15057 15144 -7511.7    15023 19.126      1  1.224e-05
## mod.herb.abun4 17 15057 15144 -7511.7    15023  0.000      0          1
## mod.herb.abun3 28 14985 15127 -7464.4    14929 94.513     11  2.164e-15
##                   
## mod.herb.abun5    
## mod.herb.abun1    
## mod.herb.abun2 ***
## mod.herb.abun4    
## mod.herb.abun3 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.herb.abun <- mod.herb.abun3
mod.herb.abun.diags <- augment(mod.herb.abun)
qplot(data=mod.herb.abun.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.herb.abun.diags, x= treatment, y=.wtres, geom="boxplot", colour=census) 

sjp.lmer(mod.herb.abun, type='fe')

sjp.lmer(mod.herb.abun, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.herb.abun, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## no.ind. ~ site + census * (treatment + exclosure) + (1 | plot/quad.unique)
##    Data: herb.abundat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##  14984.8  15127.1  -7464.4  14928.8     1162 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2136 -1.5413 -0.1873  1.3159 12.0571 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.1606   0.4007  
##  plot             (Intercept) 0.0181   0.1345  
## Number of obs: 1190, groups:  quad.unique:plot, 300; plot, 6
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.26967    0.12857   25.43  < 2e-16 ***
## site2                 -0.45906    0.14629   -3.14 0.001701 ** 
## site3                  0.02363    0.14620    0.16 0.871598    
## census16fa             0.60769    0.03466   17.53  < 2e-16 ***
## census16sp             2.10014    0.02968   70.76  < 2e-16 ***
## census16su             0.99470    0.03287   30.27  < 2e-16 ***
## treatmentF             0.19784    0.06787    2.91 0.003559 ** 
## treatmentI            -0.04490    0.07677   -0.58 0.558598    
## treatmentS            -0.87786    0.08040  -10.92  < 2e-16 ***
## treatmentW             0.03222    0.07795    0.41 0.679322    
## exclosure1            -0.10268    0.12184   -0.84 0.399373    
## census16fa:treatmentF  0.16326    0.04470    3.65 0.000260 ***
## census16sp:treatmentF  0.10428    0.03864    2.70 0.006963 ** 
## census16su:treatmentF  0.09675    0.04274    2.26 0.023608 *  
## census16fa:treatmentI  0.16519    0.04524    3.65 0.000261 ***
## census16sp:treatmentI  0.06715    0.03917    1.71 0.086438 .  
## census16su:treatmentI  0.06194    0.04312    1.44 0.150886    
## census16fa:treatmentS  0.13185    0.05613    2.35 0.018816 *  
## census16sp:treatmentS  0.24302    0.04844    5.02 5.26e-07 ***
## census16su:treatmentS  0.17013    0.05305    3.21 0.001341 ** 
## census16fa:treatmentW  0.09653    0.04497    2.15 0.031819 *  
## census16sp:treatmentW  0.05238    0.03872    1.35 0.176191    
## census16su:treatmentW  0.06326    0.04273    1.48 0.138758    
## census16fa:exclosure1  0.08300    0.03020    2.75 0.005993 ** 
## census16sp:exclosure1  0.14595    0.02622    5.57 2.61e-08 ***
## census16su:exclosure1  0.13492    0.02887    4.67 2.96e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 26 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## pesticide
herbpest.abundat <- subset(herb.abundat, treatment !='S')
herbpest.abundat$treatment <- factor(herbpest.abundat$treatment, levels=c('W', 'F', 'I', 'C'))

mod.herbpest.abun <- glmer(no.ind. ~ site + census*(treatment+exclosure)+(1|plot/quad.unique),
                           data=herbpest.abundat, family=poisson,
                           control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
mod.herbpest.abun.diags <- augment(mod.herbpest.abun)
qplot(data=mod.herbpest.abun.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

sjp.lmer(mod.herbpest.abun, type='fe')

sjp.lmer(mod.herbpest.abun, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.herbpest.abun, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## no.ind. ~ site + census * (treatment + exclosure) + (1 | plot/quad.unique)
##    Data: herbpest.abundat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##  12402.4  12519.1  -6177.2  12354.4      930 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -6.274 -1.608 -0.227  1.419 12.320 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.128522 0.35850 
##  plot             (Intercept) 0.009932 0.09966 
## Number of obs: 954, groups:  quad.unique:plot, 240; plot, 6
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.2976809  0.1060704   31.09  < 2e-16 ***
## site2                 -0.4365160  0.1150214   -3.80 0.000148 ***
## site3                 -0.0703774  0.1149445   -0.61 0.540357    
## census16fa             0.7137950  0.0354978   20.11  < 2e-16 ***
## census16sp             2.1778576  0.0307451   70.84  < 2e-16 ***
## census16su             1.0736239  0.0341081   31.48  < 2e-16 ***
## treatmentF             0.1346087  0.0729559    1.85 0.065027 .  
## treatmentI            -0.0730403  0.0749342   -0.97 0.329696    
## treatmentC            -0.0523905  0.0728981   -0.72 0.472337    
## exclosure1            -0.0430068  0.0972744   -0.44 0.658403    
## census16fa:treatmentF  0.0663157  0.0448618    1.48 0.139348    
## census16sp:treatmentF  0.0508522  0.0390387    1.30 0.192708    
## census16su:treatmentF  0.0324067  0.0431566    0.75 0.452707    
## census16fa:treatmentI  0.0690597  0.0453911    1.52 0.128150    
## census16sp:treatmentI  0.0156481  0.0395475    0.40 0.692342    
## census16su:treatmentI -0.0006187  0.0435365   -0.01 0.988661    
## census16fa:treatmentC -0.0968295  0.0449688   -2.15 0.031298 *  
## census16sp:treatmentC -0.0531242  0.0387193   -1.37 0.170053    
## census16su:treatmentC -0.0639722  0.0427254   -1.50 0.134318    
## census16fa:exclosure1  0.0622010  0.0318671    1.95 0.050952 .  
## census16sp:exclosure1  0.0936266  0.0276635    3.38 0.000713 ***
## census16su:exclosure1  0.1012586  0.0305098    3.32 0.000904 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## snow
herbsnow.abundat <- subset(herb.abundat, treatment =='C' | treatment =='S')
mod.herbsnow.abun <- glmer(no.ind. ~ site + census*treatment + exclosure+(1|plot/quad.unique),
                           data=herbsnow.abundat, family=poisson,
                           control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.herbsnow.abun.diags <- augment(mod.herbsnow.abun)
qplot(data=mod.herbsnow.abun.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

sjp.lmer(mod.herbsnow.abun, type='fe')

sjp.lmer(mod.herbsnow.abun, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

### individual species test
## without site as variables
## sample size data frame
sampsize.herbind <- as.data.frame(table(herbdat$species, herbdat$treatment, herbdat$exclosure))
colnames(sampsize.herbind) <- c('species', 'treatment', 'exclosure', 'no.quad')

sampsize.herbind <- sampsize.herbind %>% 
  spread(treatment, no.quad) 

## sample size for pesticide and snow data
sampsize.herbind.pest <- sampsize.herbind[, !(colnames(sampsize.herbind) %in% c("S"))]
sampsize.herbind.snow <- sampsize.herbind[, !(colnames(sampsize.herbind) %in% c("F",'I','W'))]

sel.sp.herbpest <- unique(sampsize.herbind.pest$species[sampsize.herbind.pest$W >= 3 & 
                        (sampsize.herbind.pest$F >= 3 | sampsize.herbind.pest$I >= 3)])
sel.sp.herbsnow <- unique(sampsize.herbind.snow$species[sampsize.herbind.snow$C >= 3 & 
                                                   sampsize.herbind.snow$S >= 3])

####################
## model with census

sampsize.herbind <- as.data.frame(table(herbdat$species, herbdat$treatment, herbdat$exclosure, herbdat$census))
colnames(sampsize.herbind) <- c('species', 'treatment', 'exclosure', 'census', 'no.quad')

sampsize.herbind <- sampsize.herbind %>% 
  spread(treatment, no.quad) 

## sample size for pesticide and snow data
sampsize.herbind.pest <- sampsize.herbind[, !(colnames(sampsize.herbind) %in% c("S"))]
sampsize.herbind.snow <- sampsize.herbind[, !(colnames(sampsize.herbind) %in% c("F",'I','W'))]

sel.sp.herbpest <- unique(sampsize.herbind.pest$species[sampsize.herbind.pest$W >= 3 & 
                        (sampsize.herbind.pest$F >= 3 | sampsize.herbind.pest$I >= 3)])
sel.sp.herbsnow <- unique(sampsize.herbind.snow$species[sampsize.herbind.snow$C >= 3 & 
                                                   sampsize.herbind.snow$S >= 3])

indivmods.herbpest <- sapply(sel.sp.herbpest, function(sp){
  print(sp)
  spdat <- filter(herbpest, species == sp)
  spdat <- droplevels(spdat)
  # random effects: remove quad.unique
  if (length(unique(spdat$exclosure)) == 2 & length(unique(spdat$census)) >= 2){
    mod <- glmer(no.ind. ~ census + treatment + exclosure + (1|plot/quad.unique),
                 data=spdat,family=poisson,
                 control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
  } else if (length(unique(spdat$exclosure)) < 2 & length(unique(spdat$census)) >=2) {
    mod <- glmer(no.ind. ~ census + treatment + (1|plot/quad.unique),
                 data=spdat,family=poisson,
                 control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
  } else if (length(unique(spdat$exclosure)) == 2 & length(unique(spdat$census)) < 2) {
    mod <- glmer(no.ind. ~ treatment + exclosure + (1|plot/quad.unique),
                 data=spdat, family=poisson,
                 control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
  } else {
    mod <- glmer(no.ind. ~ treatment + (1|plot/quad.unique),
                 data=spdat, family=poisson,
                 control=glmerControl(optimizer="optimx", 
                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
  }
  
}, simplify=FALSE)
## [1] Adenophora remotiflora
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Adonis amurensis
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Adoxa moschatellina
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0508559 (tol =
## 0.001, component 1)
## [1] Allium monanthum
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Anemone amurensis
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Anemone baicalensis
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Anemone raddeana
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Anemone spp.
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Arisaema amurense
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Athyrium brevifrons
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Brachybotrys paridiformis
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Cardamine leucantha
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] Carex pilosa
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Carex remotiuscula
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] Circaea cordata
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Circaea lutetiana
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Corydalis repens
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Corydalis turtschaninovii
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Diarrhena manshurica
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Enemion raddeanum
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Equisetum hyemale
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Eranthis stellata
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Filipendula palmata
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Gagea pauciflora
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Gagea triflora
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Galium triflorum
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Impatiens noli-tangere
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Lamium barbatum
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Lilium distichum
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] Maianthemum bifolium
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Meehania urticifolia
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Osmorhiza aristata
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Oxalis acetosella subsp. Griffithii
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Parasenecio hastatus
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] Paris verticillata
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0436602 (tol =
## 0.001, component 1)
## [1] Phryma leptostachya subsp. asiatica
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] Pseudostellaria sylvatica
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Ranunculus japonicus
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Saussurea manshurica
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Smilacina japonica
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Thalictrum spp.
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] unknown spp2
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] unknown spp4
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
## [1] Viola verecunda
## 82 Levels: Aconitum kusnezoffii Aconitum sczukinii ... Viola verecunda
### plot residuals
lapply(indivmods.herbpest, function(mod){
  plot(mod, type=c('p', 'smooth'))})
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

## 
## [[33]]

## 
## [[34]]

## 
## [[35]]

## 
## [[36]]

## 
## [[37]]

## 
## [[38]]

## 
## [[39]]

## 
## [[40]]

## 
## [[41]]

## 
## [[42]]

## 
## [[43]]

## 
## [[44]]

lapply(indivmods.herbpest, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## [[1]]
## [[1]]$`quad.unique:plot`

## 
## [[1]]$plot

## 
## 
## [[2]]
## [[2]]$`quad.unique:plot`

## 
## [[2]]$plot

## 
## 
## [[3]]
## [[3]]$`quad.unique:plot`

## 
## [[3]]$plot

## 
## 
## [[4]]
## [[4]]$`quad.unique:plot`

## 
## [[4]]$plot

## 
## 
## [[5]]
## [[5]]$`quad.unique:plot`

## 
## [[5]]$plot

## 
## 
## [[6]]
## [[6]]$`quad.unique:plot`

## 
## [[6]]$plot

## 
## 
## [[7]]
## [[7]]$`quad.unique:plot`

## 
## [[7]]$plot

## 
## 
## [[8]]
## [[8]]$`quad.unique:plot`

## 
## [[8]]$plot

## 
## 
## [[9]]
## [[9]]$`quad.unique:plot`

## 
## [[9]]$plot

## 
## 
## [[10]]
## [[10]]$`quad.unique:plot`

## 
## [[10]]$plot

## 
## 
## [[11]]
## [[11]]$`quad.unique:plot`

## 
## [[11]]$plot

## 
## 
## [[12]]
## [[12]]$`quad.unique:plot`

## 
## [[12]]$plot

## 
## 
## [[13]]
## [[13]]$`quad.unique:plot`

## 
## [[13]]$plot

## 
## 
## [[14]]
## [[14]]$`quad.unique:plot`

## 
## [[14]]$plot

## 
## 
## [[15]]
## [[15]]$`quad.unique:plot`

## 
## [[15]]$plot

## 
## 
## [[16]]
## [[16]]$`quad.unique:plot`

## 
## [[16]]$plot

## 
## 
## [[17]]
## [[17]]$`quad.unique:plot`

## 
## [[17]]$plot

## 
## 
## [[18]]
## [[18]]$`quad.unique:plot`

## 
## [[18]]$plot

## 
## 
## [[19]]
## [[19]]$`quad.unique:plot`

## 
## [[19]]$plot

## 
## 
## [[20]]
## [[20]]$`quad.unique:plot`

## 
## [[20]]$plot

## 
## 
## [[21]]
## [[21]]$`quad.unique:plot`

## 
## [[21]]$plot

## 
## 
## [[22]]
## [[22]]$`quad.unique:plot`

## 
## [[22]]$plot

## 
## 
## [[23]]
## [[23]]$`quad.unique:plot`

## 
## [[23]]$plot

## 
## 
## [[24]]
## [[24]]$`quad.unique:plot`

## 
## [[24]]$plot

## 
## 
## [[25]]
## [[25]]$`quad.unique:plot`

## 
## [[25]]$plot

## 
## 
## [[26]]
## [[26]]$`quad.unique:plot`

## 
## [[26]]$plot

## 
## 
## [[27]]
## [[27]]$`quad.unique:plot`

## 
## [[27]]$plot

## 
## 
## [[28]]
## [[28]]$`quad.unique:plot`

## 
## [[28]]$plot

## 
## 
## [[29]]
## [[29]]$`quad.unique:plot`

## 
## [[29]]$plot

## 
## 
## [[30]]
## [[30]]$`quad.unique:plot`

## 
## [[30]]$plot

## 
## 
## [[31]]
## [[31]]$`quad.unique:plot`

## 
## [[31]]$plot

## 
## 
## [[32]]
## [[32]]$`quad.unique:plot`

## 
## [[32]]$plot

## 
## 
## [[33]]
## [[33]]$`quad.unique:plot`

## 
## [[33]]$plot

## 
## 
## [[34]]
## [[34]]$`quad.unique:plot`

## 
## [[34]]$plot

## 
## 
## [[35]]
## [[35]]$`quad.unique:plot`

## 
## [[35]]$plot

## 
## 
## [[36]]
## [[36]]$`quad.unique:plot`

## 
## [[36]]$plot

## 
## 
## [[37]]
## [[37]]$`quad.unique:plot`

## 
## [[37]]$plot

## 
## 
## [[38]]
## [[38]]$`quad.unique:plot`

## 
## [[38]]$plot

## 
## 
## [[39]]
## [[39]]$`quad.unique:plot`

## 
## [[39]]$plot

## 
## 
## [[40]]
## [[40]]$`quad.unique:plot`

## 
## [[40]]$plot

## 
## 
## [[41]]
## [[41]]$`quad.unique:plot`

## 
## [[41]]$plot

## 
## 
## [[42]]
## [[42]]$`quad.unique:plot`

## 
## [[42]]$plot

## 
## 
## [[43]]
## [[43]]$`quad.unique:plot`

## 
## [[43]]$plot

## 
## 
## [[44]]
## [[44]]$`quad.unique:plot`

## 
## [[44]]$plot

lapply(indivmods.herbpest, function(x) summary(x, correlation=FALSE))
## [[1]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    396.3    418.5   -188.2    376.3       58 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3896 -0.7003 -0.3194  0.5528  2.1206 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.562    0.7497  
##  plot             (Intercept) 0.000    0.0000  
## Number of obs: 68, groups:  quad.unique:plot, 34; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.1932     0.7948  -0.243 0.807894    
## census16fa    1.4054     0.7264   1.935 0.053011 .  
## census16sp    2.2337     0.7194   3.105 0.001901 ** 
## census16su    2.6442     0.7180   3.683 0.000231 ***
## treatmentF    0.3359     0.3828   0.878 0.380196    
## treatmentI   -0.2622     0.4426  -0.592 0.553666    
## treatmentC    0.2975     0.4219   0.705 0.480689    
## exclosure1   -1.1350     0.3280  -3.461 0.000539 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[2]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    231.2    243.9   -108.6    217.2       38 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.24302 -0.68936  0.03475  0.51888  1.15953 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2019   0.4493  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 45, groups:  quad.unique:plot, 45; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.06652    0.24667   4.324 1.53e-05 ***
## treatmentF   0.59801    0.29169   2.050   0.0403 *  
## treatmentI   0.54540    0.30607   1.782   0.0748 .  
## treatmentC   0.39864    0.28265   1.410   0.1584    
## exclosure1  -0.03868    0.19983  -0.194   0.8465    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[3]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1453.7   1487.6   -717.8   1435.7      311 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61687 -0.53196 -0.08948  0.37939  2.68756 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.246061 0.49605 
##  plot             (Intercept) 0.002631 0.05129 
## Number of obs: 320, groups:  quad.unique:plot, 163; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.051998   0.001418    36.7   <2e-16 ***
## census16sp  1.214466   0.001418   856.5   <2e-16 ***
## census16su  1.000188   0.001418   705.5   <2e-16 ***
## treatmentF  0.251437   0.001417   177.4   <2e-16 ***
## treatmentI  0.318203   0.001417   224.5   <2e-16 ***
## treatmentC  0.154555   0.001417   109.0   <2e-16 ***
## exclosure1  0.029006   0.001418    20.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0508559 (tol = 0.001, component 1)
## 
## 
## [[4]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1391.8   1414.2   -688.9   1377.8      175 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.26042 -0.33517  0.00582  0.18330  0.33528 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.62654  0.7915  
##  plot             (Intercept) 0.06826  0.2613  
## Number of obs: 182, groups:  quad.unique:plot, 182; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.46699    0.21105  11.689   <2e-16 ***
## treatmentF   0.28454    0.18147   1.568    0.117    
## treatmentI  -0.06145    0.18214  -0.337    0.736    
## treatmentC  -0.01336    0.17877  -0.075    0.940    
## exclosure1  -0.07854    0.24936  -0.315    0.753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[5]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   2435.7   2459.9  -1210.9   2421.7      226 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75531 -0.14240  0.01566  0.08460  0.20769 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.7202   0.8487  
##  plot             (Intercept) 0.1230   0.3506  
## Number of obs: 233, groups:  quad.unique:plot, 233; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.0827     0.2394  17.053   <2e-16 ***
## treatmentF   -0.1200     0.1610  -0.745    0.456    
## treatmentI   -0.1570     0.1609  -0.976    0.329    
## treatmentC   -0.0646     0.1600  -0.404    0.686    
## exclosure1   -0.2479     0.3080  -0.805    0.421    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[6]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    791.9    817.0   -385.9    771.9       81 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8030 -1.1121 -0.3899  0.6558  5.1412 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.7416   0.8612  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 91, groups:  quad.unique:plot, 34; plot, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.4328     0.4139   1.046    0.296    
## census16fa    0.9866     0.2158   4.571 4.85e-06 ***
## census16sp    2.8923     0.1957  14.779  < 2e-16 ***
## census16su    3.2737     0.1949  16.798  < 2e-16 ***
## treatmentF    0.2980     0.4436   0.672    0.502    
## treatmentI   -0.3831     0.4168  -0.919    0.358    
## treatmentC   -0.3491     0.4157  -0.840    0.401    
## exclosure1   -0.1249     0.3181  -0.393    0.694    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[7]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    590.3    609.8   -287.2    574.3       76 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.24459 -0.28497 -0.02441  0.22054  0.45976 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.5669   0.7529  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 84, groups:  quad.unique:plot, 83; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.64643    1.07580   0.601    0.548
## census16sp   1.49300    1.05422   1.416    0.157
## treatmentF   0.13113    0.27969   0.469    0.639
## treatmentI   0.30633    0.25198   1.216    0.224
## treatmentC  -0.04721    0.25262  -0.187    0.852
## exclosure1  -0.02352    0.19811  -0.119    0.905
## 
## [[8]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    170.3    179.1    -78.1    156.3       19 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.18169 -0.40673  0.09073  0.28263  0.64728 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.4823   0.6945  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 26, groups:  quad.unique:plot, 26; plot, 2
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.61705    0.32227   5.018 5.23e-07 ***
## treatmentF   0.50351    0.42561   1.183    0.237    
## treatmentI   0.50689    0.42585   1.190    0.234    
## treatmentC  -0.43236    0.48357  -0.894    0.371    
## exclosure1   0.04706    0.33168   0.142    0.887    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[9]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    297.7    318.8   -138.9    277.7       51 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.47053 -0.51244 -0.09982  0.41978  2.07621 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.21929  0.4683  
##  plot             (Intercept) 0.09572  0.3094  
## Number of obs: 61, groups:  quad.unique:plot, 29; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.87853    0.35001   5.367 8.00e-08 ***
## census16fa  -0.97916    0.22884  -4.279 1.88e-05 ***
## census16sp  -0.52181    0.16400  -3.182  0.00146 ** 
## census16su   0.22807    0.13591   1.678  0.09333 .  
## treatmentF  -0.86181    0.30199  -2.854  0.00432 ** 
## treatmentI  -0.83382    0.45314  -1.840  0.06575 .  
## treatmentC  -1.14623    0.38258  -2.996  0.00274 ** 
## exclosure1   0.08685    0.36775   0.236  0.81331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[10]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    255.3    281.5   -117.7    235.3       91 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.53991 -0.25004 -0.19831  0.02704  1.43182 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0        0       
##  plot             (Intercept) 0        0       
## Number of obs: 101, groups:  quad.unique:plot, 50; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.15353    0.24826   0.618    0.536
## census16fa  -0.14567    0.26787  -0.544    0.587
## census16sp   0.09562    0.25385   0.377    0.706
## census16su   0.09586    0.24574   0.390    0.696
## treatmentF   0.01306    0.30137   0.043    0.965
## treatmentI   0.20481    0.25113   0.816    0.415
## treatmentC  -0.03490    0.26131  -0.134    0.894
## exclosure1   0.07936    0.20317   0.391    0.696
## 
## [[11]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   2306.3   2347.3  -1143.1   2286.3      438 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1596 -0.5486 -0.1378  0.4045  3.3871 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.6707   0.8190  
##  plot             (Intercept) 0.1958   0.4425  
## Number of obs: 448, groups:  quad.unique:plot, 157; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.56591    0.30435   1.859    0.063 .  
## census16fa   0.74335    0.07709   9.642   <2e-16 ***
## census16sp   1.28434    0.07287  17.626   <2e-16 ***
## census16su   1.27029    0.07409  17.145   <2e-16 ***
## treatmentF  -0.10269    0.20494  -0.501    0.616    
## treatmentI   0.08965    0.20174   0.444    0.657    
## treatmentC  -0.08390    0.19493  -0.430    0.667    
## exclosure1  -0.42943    0.39247  -1.094    0.274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[12]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5819.7   5867.8  -2899.8   5799.7      897 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4861 -0.8254 -0.1137  0.6283  8.0141 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2081   0.4561  
##  plot             (Intercept) 0.1155   0.3398  
## Number of obs: 907, groups:  quad.unique:plot, 237; plot, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.516592   0.209980   7.223  5.1e-13 ***
## census16fa   0.652098   0.029934  21.785  < 2e-16 ***
## census16sp   0.706036   0.029585  23.865  < 2e-16 ***
## census16su   0.529590   0.030727  17.235  < 2e-16 ***
## treatmentF   0.186653   0.088761   2.103   0.0355 *  
## treatmentI   0.175275   0.089792   1.952   0.0509 .  
## treatmentC  -0.000753   0.089653  -0.008   0.9933    
## exclosure1   0.355648   0.284645   1.249   0.2115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## [[13]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1688.4   1726.6   -834.2   1668.4      327 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5181 -0.5681 -0.2027  0.3797  3.8899 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.6483   0.8052  
##  plot             (Intercept) 0.0745   0.2729  
## Number of obs: 337, groups:  quad.unique:plot, 99; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.06121    0.24995   4.246 2.18e-05 ***
## census16fa   0.25452    0.06763   3.763 0.000168 ***
## census16sp   0.39441    0.06630   5.949 2.70e-09 ***
## census16su   0.48459    0.06543   7.406 1.30e-13 ***
## treatmentF   0.16925    0.24509   0.691 0.489839    
## treatmentI  -0.04732    0.24481  -0.193 0.846716    
## treatmentC  -0.31457    0.25351  -1.241 0.214667    
## exclosure1   0.02077    0.28664   0.072 0.942243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[14]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    457.8    484.0   -218.9    437.8       92 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1758 -0.4776 -0.2514  0.3632  3.5140 
## 
## Random effects:
##  Groups           Name        Variance  Std.Dev. 
##  quad.unique:plot (Intercept) 5.671e-01 7.531e-01
##  plot             (Intercept) 6.461e-16 2.542e-08
## Number of obs: 102, groups:  quad.unique:plot, 47; plot, 6
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.0006436  0.3817867  -0.002 0.998655    
## census16fa   0.4229696  0.1761388   2.401 0.016335 *  
## census16sp   0.6425333  0.1664055   3.861 0.000113 ***
## census16su   0.5212396  0.1760831   2.960 0.003074 ** 
## treatmentF   0.7001961  0.3739002   1.873 0.061112 .  
## treatmentI   0.4659222  0.4066670   1.146 0.251915    
## treatmentC   0.2652807  0.4591589   0.578 0.563431    
## exclosure1  -0.1844836  0.2759277  -0.669 0.503755    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## [[15]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    516.8    540.4   -249.4    498.8       93 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4146 -0.4890 -0.2041  0.3228  3.3416 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.5469   0.7395  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 102, groups:  quad.unique:plot, 85; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.24188    1.08625  -0.223   0.8238  
## census16fa   1.26387    1.08452   1.165   0.2439  
## census16su   1.54344    1.07608   1.434   0.1515  
## treatmentF   0.11730    0.27892   0.420   0.6741  
## treatmentI  -0.04357    0.28990  -0.150   0.8805  
## treatmentC  -0.19707    0.28332  -0.696   0.4867  
## exclosure1  -0.37475    0.20785  -1.803   0.0714 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[16]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    448.0    474.3   -214.0    428.0       92 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5822 -0.5450 -0.1479  0.3632  1.8044 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.3729   0.6107  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 102, groups:  quad.unique:plot, 57; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.8319     0.2383   3.492  0.00048 ***
## census16fa   -0.4416     0.2422  -1.823  0.06829 .  
## census16sp   -1.1130     0.2607  -4.270 1.95e-05 ***
## census16su    0.2254     0.1224   1.841  0.06564 .  
## treatmentF    0.4427     0.3145   1.408  0.15923    
## treatmentI    0.2382     0.2774   0.859  0.39055    
## treatmentC    0.1388     0.3099   0.448  0.65424    
## exclosure1   -0.1934     0.2128  -0.909  0.36346    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[17]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1871.1   1894.7   -928.5   1857.1      209 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.82777 -0.23149  0.02501  0.15341  0.41889 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.4487   0.6699  
##  plot             (Intercept) 0.2413   0.4912  
## Number of obs: 216, groups:  quad.unique:plot, 216; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.9102     0.3058   9.518   <2e-16 ***
## treatmentF    0.2054     0.1387   1.481    0.139    
## treatmentI    0.1852     0.1374   1.347    0.178    
## treatmentC    0.2019     0.1380   1.463    0.144    
## exclosure1    0.1219     0.4129   0.295    0.768    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[18]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1147.8   1173.9   -565.9   1131.8      185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.35606 -0.44858 -0.09649  0.27186  1.01838 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.4113   0.6413  
##  plot             (Intercept) 0.1701   0.4125  
## Number of obs: 193, groups:  quad.unique:plot, 189; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 0.496709   0.438669   1.132  0.25750   
## census16sp  1.005545   0.339429   2.962  0.00305 **
## treatmentF  0.193078   0.157676   1.224  0.22075   
## treatmentI  0.004239   0.160718   0.026  0.97896   
## treatmentC  0.316229   0.159704   1.980  0.04769 * 
## exclosure1  0.120181   0.361494   0.332  0.73954   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[19]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    420.1    446.8   -200.0    400.1       97 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2358 -0.5267 -0.1009  0.4238  2.1813 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.09759  0.3124  
##  plot             (Intercept) 0.07865  0.2804  
## Number of obs: 107, groups:  quad.unique:plot, 38; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.94602    0.29890   3.165 0.001551 ** 
## census16fa   0.43601    0.16589   2.628 0.008582 ** 
## census16sp   0.34518    0.16475   2.095 0.036156 *  
## census16su   0.57385    0.15736   3.647 0.000266 ***
## treatmentF  -0.78230    0.26096  -2.998 0.002719 ** 
## treatmentI   0.10917    0.21307   0.512 0.608388    
## treatmentC  -0.58940    0.25678  -2.295 0.021714 *  
## exclosure1  -0.09217    0.30322  -0.304 0.761151    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[20]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    398.0    420.1   -190.0    380.0       77 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3837 -0.5763 -0.2547  0.3128  1.8627 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.4198   0.6479  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 86, groups:  quad.unique:plot, 47; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.8803     0.3536   2.490  0.01279 * 
## census16sp    0.7994     0.2541   3.145  0.00166 **
## census16su    0.7488     0.2519   2.973  0.00295 **
## treatmentF   -0.7385     0.3593  -2.056  0.03982 * 
## treatmentI   -0.7533     0.3169  -2.377  0.01745 * 
## treatmentC   -0.7411     0.3341  -2.218  0.02656 * 
## exclosure1    0.0936     0.2492   0.376  0.70716   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[21]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    812.3    847.7   -396.2    792.3      243 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3643 -0.4603 -0.1474  0.3077  4.4097 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.1      0.3163  
##  plot             (Intercept) 0.0      0.0000  
## Number of obs: 253, groups:  quad.unique:plot, 106; plot, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.406547   0.168164   2.418  0.01563 * 
## census16fa   0.426732   0.145646   2.930  0.00339 **
## census16sp   0.045427   0.168414   0.270  0.78736   
## census16su   0.347470   0.147283   2.359  0.01831 * 
## treatmentF  -0.165827   0.175827  -0.943  0.34562   
## treatmentI  -0.006441   0.151502  -0.042  0.96609   
## treatmentC   0.158244   0.155209   1.020  0.30794   
## exclosure1  -0.239882   0.121035  -1.982  0.04749 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[22]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1092.9   1118.0   -538.4   1076.9      163 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.36956 -0.44994 -0.02384  0.25681  2.23279 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.58779  0.7667  
##  plot             (Intercept) 0.02272  0.1507  
## Number of obs: 171, groups:  quad.unique:plot, 158; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.4316     0.1958  12.421  < 2e-16 ***
## census16su   -1.0351     0.1590  -6.509 7.56e-11 ***
## treatmentF   -0.4905     0.1995  -2.458   0.0140 *  
## treatmentI   -0.3097     0.1954  -1.585   0.1130    
## treatmentC   -0.1698     0.1984  -0.856   0.3921    
## exclosure1   -0.3483     0.1947  -1.789   0.0737 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[23]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1116.2   1149.4   -548.1   1096.2      195 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1002 -0.7797 -0.2281  0.5473  5.9680 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.35126  0.5927  
##  plot             (Intercept) 0.04874  0.2208  
## Number of obs: 205, groups:  quad.unique:plot, 63; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.91113    0.25093   3.631 0.000282 ***
## census16fa   0.25810    0.09705   2.660 0.007825 ** 
## census16sp   0.91129    0.08503  10.717  < 2e-16 ***
## census16su   0.19852    0.09572   2.074 0.038075 *  
## treatmentF   0.37161    0.23466   1.584 0.113282    
## treatmentI   0.26893    0.23007   1.169 0.242440    
## treatmentC   0.23311    0.24399   0.955 0.339372    
## exclosure1  -0.26008    0.28984  -0.897 0.369550    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[24]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    677.3    694.4   -331.6    663.3       79 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.83772 -0.31638 -0.00503  0.13073  0.28653 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.9296   0.9642  
##  plot             (Intercept) 0.1771   0.4209  
## Number of obs: 86, groups:  quad.unique:plot, 86; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.47709    0.40236   6.156 7.44e-10 ***
## treatmentF   0.12513    0.32110   0.390    0.697    
## treatmentI   0.00742    0.32446   0.023    0.982    
## treatmentC  -0.09743    0.31091  -0.313    0.754    
## exclosure1  -0.48860    0.45873  -1.065    0.287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[25]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    917.8    935.9   -451.9    903.8       90 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.02578 -0.22425  0.02181  0.08840  0.22665 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.9977   0.9989  
##  plot             (Intercept) 0.2590   0.5089  
## Number of obs: 97, groups:  quad.unique:plot, 97; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.37986    0.43852   7.707 1.28e-14 ***
## treatmentF  -0.03696    0.30356  -0.122    0.903    
## treatmentI   0.01052    0.30494   0.034    0.972    
## treatmentC  -0.13928    0.29507  -0.472    0.637    
## exclosure1  -0.63895    0.51862  -1.232    0.218    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[26]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1306.1   1341.1   -643.1   1286.1      235 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5698 -0.6979 -0.2716  0.5267  3.6836 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.5044   0.7102  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 245, groups:  quad.unique:plot, 73; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.86255    0.21419   4.027 5.65e-05 ***
## census16fa   0.57119    0.08653   6.601 4.08e-11 ***
## census16sp   1.01998    0.08025  12.709  < 2e-16 ***
## census16su   0.66204    0.08635   7.667 1.76e-14 ***
## treatmentF   0.25668    0.25102   1.023    0.307    
## treatmentI   0.30859    0.25626   1.204    0.229    
## treatmentC   0.14443    0.25786   0.560    0.575    
## exclosure1  -0.24376    0.18144  -1.343    0.179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[27]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1549.2   1586.3   -764.6   1529.2      292 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1426 -0.6623 -0.2713  0.5493  4.1255 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2614   0.5112  
##  plot             (Intercept) 0.1027   0.3205  
## Number of obs: 302, groups:  quad.unique:plot, 115; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.73492    0.31080   2.365    0.018 *  
## census16fa   0.21424    0.15315   1.399    0.162    
## census16sp   0.81354    0.14315   5.683 1.32e-08 ***
## census16su   0.91411    0.14198   6.438 1.21e-10 ***
## treatmentF  -0.06830    0.15528  -0.440    0.660    
## treatmentI  -0.21628    0.15811  -1.368    0.171    
## treatmentC  -0.21779    0.15106  -1.442    0.149    
## exclosure1   0.06945    0.34084   0.204    0.839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[28]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    252.2    271.5   -116.1    232.2       41 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.09730 -0.55654 -0.08424  0.32248  1.64176 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.3401   0.5831  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 51, groups:  quad.unique:plot, 35; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.72850    0.40342   4.285 1.83e-05 ***
## census16fa  -0.46234    0.34695  -1.333   0.1827    
## census16sp   0.36255    0.60775   0.597   0.5508    
## census16su   0.01212    0.29199   0.042   0.9669    
## treatmentF  -0.37011    0.35388  -1.046   0.2956    
## treatmentI  -0.06169    0.45361  -0.136   0.8918    
## treatmentC   0.07508    0.34541   0.217   0.8279    
## exclosure1  -0.54254    0.26530  -2.045   0.0409 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[29]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    466.6    490.0   -224.3    448.6       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6922 -0.5343 -0.3512  0.3043  2.4551 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.5122   0.7157  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 99, groups:  quad.unique:plot, 70; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.32745    0.41462  -0.790    0.430    
## census16sp   1.21991    0.30538   3.995 6.48e-05 ***
## census16su   1.26755    0.30545   4.150 3.33e-05 ***
## treatmentF  -0.18670    0.32422  -0.576    0.565    
## treatmentI  -0.23109    0.33265  -0.695    0.487    
## treatmentC  -0.07917    0.31578  -0.251    0.802    
## exclosure1   0.12679    0.22729   0.558    0.577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## [[30]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    943.6    976.8   -461.8    923.6      193 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5080 -0.5334 -0.1544  0.2482  2.3298 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.6794   0.8243  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 203, groups:  quad.unique:plot, 93; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.19556    0.74334  -1.608 0.107758    
## census16fa   1.78384    0.71691   2.488 0.012838 *  
## census16sp   2.45628    0.71479   3.436 0.000590 ***
## census16su   2.37541    0.71492   3.323 0.000892 ***
## treatmentF  -0.35364    0.27597  -1.281 0.200040    
## treatmentI  -0.19453    0.27471  -0.708 0.478870    
## treatmentC   0.07528    0.25437   0.296 0.767278    
## exclosure1   0.06576    0.19683   0.334 0.738322    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[31]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5094.6   5140.8  -2537.3   5074.6      736 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8240 -0.7520 -0.1286  0.6838  5.6426 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.8039   0.8966  
##  plot             (Intercept) 0.1097   0.3312  
## Number of obs: 746, groups:  quad.unique:plot, 218; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.59569    0.24102   6.621 3.57e-11 ***
## census16fa   0.67614    0.03431  19.710  < 2e-16 ***
## census16sp   1.01412    0.03288  30.846  < 2e-16 ***
## census16su   0.76758    0.03441  22.306  < 2e-16 ***
## treatmentF  -0.15405    0.17390  -0.886    0.376    
## treatmentI  -0.25253    0.17854  -1.414    0.157    
## treatmentC  -0.07753    0.16991  -0.456    0.648    
## exclosure1  -0.02563    0.29868  -0.086    0.932    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[32]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   2622.7   2664.5  -1301.3   2602.7      473 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5130 -0.7168 -0.1679  0.4744  4.5669 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.45287  0.673   
##  plot             (Intercept) 0.02857  0.169   
## Number of obs: 483, groups:  quad.unique:plot, 173; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.67036    0.16838   3.981 6.86e-05 ***
## census16fa   1.01894    0.07771  13.112  < 2e-16 ***
## census16sp   0.80273    0.08001  10.033  < 2e-16 ***
## census16su   1.17841    0.07665  15.375  < 2e-16 ***
## treatmentF   0.03023    0.15825   0.191    0.848    
## treatmentI  -0.12507    0.16025  -0.780    0.435    
## treatmentC  -0.04981    0.15411  -0.323    0.747    
## exclosure1  -0.08124    0.18054  -0.450    0.653    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[33]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   2576.2   2615.5  -1278.1   2556.2      364 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1081 -0.8340 -0.2325  0.5387  7.6897 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.7443   0.8627  
##  plot             (Intercept) 0.1065   0.3264  
## Number of obs: 374, groups:  quad.unique:plot, 121; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.80501    0.25732   7.015  2.3e-12 ***
## census16fa   0.52834    0.04380  12.063  < 2e-16 ***
## census16sp  -0.59608    0.05910 -10.087  < 2e-16 ***
## census16su   0.01491    0.04904   0.304   0.7611    
## treatmentF  -0.52044    0.24540  -2.121   0.0339 *  
## treatmentI  -0.25108    0.23535  -1.067   0.2861    
## treatmentC  -0.00968    0.22498  -0.043   0.9657    
## exclosure1  -0.21234    0.34632  -0.613   0.5398    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[34]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    246.9    267.5   -113.5    226.9       48 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1712 -0.4479 -0.1177  0.2456  1.8952 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2385   0.4883  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 58, groups:  quad.unique:plot, 28; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.97986    0.58991  -1.661 0.096709 .  
## census16fa   1.42909    0.53144   2.689 0.007165 ** 
## census16sp   1.58312    0.52823   2.997 0.002726 ** 
## census16su   1.91969    0.52285   3.672 0.000241 ***
## treatmentF   0.75901    0.37978   1.999 0.045657 *  
## treatmentI   0.05891    0.39017   0.151 0.879984    
## treatmentC   0.21187    0.33084   0.640 0.521909    
## exclosure1   0.18100    0.27437   0.660 0.509450    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## [[35]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    769.6    799.6   -374.8    749.6      138 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87801 -0.59860 -0.02626  0.31673  1.88947 
## 
## Random effects:
##  Groups           Name        Variance  Std.Dev. 
##  quad.unique:plot (Intercept) 4.557e-01 0.6750686
##  plot             (Intercept) 2.033e-08 0.0001426
## Number of obs: 148, groups:  quad.unique:plot, 106; plot, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.039266   0.001637    24.0   <2e-16 ***
## census16fa  -0.037781   0.001635   -23.1   <2e-16 ***
## census16sp   1.606745   0.001637   981.5   <2e-16 ***
## census16su   0.290372   0.001635   177.6   <2e-16 ***
## treatmentF  -0.014153   0.001636    -8.7   <2e-16 ***
## treatmentI   0.183950   0.001636   112.4   <2e-16 ***
## treatmentC   0.038881   0.001636    23.8   <2e-16 ***
## exclosure1   0.003387   0.001637     2.1   0.0385 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0436602 (tol = 0.001, component 1)
## 
## 
## [[36]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    215.8    233.2    -98.9    197.8       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9748 -0.5488 -0.2258  0.2114  1.5869 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.1921   0.4383  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 51, groups:  quad.unique:plot, 41; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.72428    0.49952   1.450   0.1471  
## census16fa   0.58484    0.42723   1.369   0.1710  
## census16su   0.34953    0.43729   0.799   0.4241  
## treatmentF   0.10343    0.40117   0.258   0.7965  
## treatmentI  -0.06969    0.33273  -0.210   0.8341  
## treatmentC  -0.03335    0.34747  -0.096   0.9235  
## exclosure1  -0.55898    0.23687  -2.360   0.0183 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## [[37]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1443.7   1482.0   -711.8   1423.7      330 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9036 -0.5516 -0.1287  0.3094  3.3025 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2945   0.5426  
##  plot             (Intercept) 0.1098   0.3314  
## Number of obs: 340, groups:  quad.unique:plot, 170; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.20104    1.07408   0.187    0.852
## census16fa   0.17542    1.06211   0.165    0.869
## census16sp   0.85556    1.05832   0.808    0.419
## census16su   0.97918    1.05906   0.925    0.355
## treatmentF  -0.02294    0.15363  -0.149    0.881
## treatmentI  -0.18000    0.14700  -1.224    0.221
## treatmentC  -0.04480    0.14935  -0.300    0.764
## exclosure1   0.03831    0.29237   0.131    0.896
## 
## [[38]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    513.1    534.3   -248.6    497.1       96 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1356 -0.4775 -0.1764  0.3991  1.1174 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.28393  0.5329  
##  plot             (Intercept) 0.04642  0.2155  
## Number of obs: 104, groups:  quad.unique:plot, 103; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.49457    0.50856   0.972    0.331
## census16sp   0.67473    0.46207   1.460    0.144
## treatmentF   0.06742    0.22538   0.299    0.765
## treatmentI   0.14440    0.20139   0.717    0.473
## treatmentC   0.14406    0.20546   0.701    0.483
## exclosure1   0.12906    0.23017   0.561    0.575
## 
## [[39]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    210.3    231.1    -95.2    190.3       49 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8695 -0.4751 -0.1201  0.2784  2.0151 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.1296   0.3601  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 59, groups:  quad.unique:plot, 26; plot, 6
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.2481969  0.6873264  -0.361   0.7180  
## census16fa   0.4961348  0.6463335   0.768   0.4427  
## census16sp   1.1263778  0.6302712   1.787   0.0739 .
## census16su   1.0602241  0.6425455   1.650   0.0989 .
## treatmentF  -0.0005728  0.3284512  -0.002   0.9986  
## treatmentI  -0.4627426  0.4615937  -1.002   0.3161  
## treatmentC   0.2566551  0.3293613   0.779   0.4358  
## exclosure1  -0.0647791  0.2408050  -0.269   0.7879  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[40]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    284.2    304.1   -132.1    264.2       44 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.91485 -0.46320 -0.04718  0.39584  1.29404 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 1.02560  1.0127  
##  plot             (Intercept) 0.01845  0.1358  
## Number of obs: 54, groups:  quad.unique:plot, 27; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.3287     1.1320  -1.174  0.24046   
## census16fa    1.5245     1.0486   1.454  0.14601   
## census16sp    3.2876     1.0183   3.229  0.00124 **
## census16su    2.9149     1.0295   2.831  0.00464 **
## treatmentF    0.1851     0.6727   0.275  0.78326   
## treatmentI   -0.9168     1.2068  -0.760  0.44744   
## treatmentC   -0.4440     0.7843  -0.566  0.57133   
## exclosure1   -0.6273     0.5622  -1.116  0.26456   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[41]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    326.8    347.8   -154.4    308.8       67 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0176 -0.4514 -0.1487  0.2348  1.6057 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.40685  0.6379  
##  plot             (Intercept) 0.06955  0.2637  
## Number of obs: 76, groups:  quad.unique:plot, 38; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.6372     0.3908   1.631   0.1030  
## census16sp    0.4353     0.1953   2.229   0.0258 *
## census16su    0.4842     0.1955   2.477   0.0132 *
## treatmentF    0.1058     0.3941   0.269   0.7883  
## treatmentI    0.3676     0.3821   0.962   0.3359  
## treatmentC   -0.3545     0.3757  -0.944   0.3453  
## exclosure1   -0.1642     0.3819  -0.430   0.6672  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[42]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    302.1    326.2   -141.1    282.1       72 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.93573 -0.46663 -0.07868  0.33451  1.78963 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.2062   0.4541  
##  plot             (Intercept) 0.0000   0.0000  
## Number of obs: 82, groups:  quad.unique:plot, 46; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.03789    0.49821   0.076   0.9394  
## census16fa   0.36163    0.43452   0.832   0.4053  
## census16sp   0.30640    0.43271   0.708   0.4789  
## census16su   0.53355    0.42953   1.242   0.2142  
## treatmentF   0.33052    0.34681   0.953   0.3406  
## treatmentI   0.48659    0.30043   1.620   0.1053  
## treatmentC   0.53673    0.31669   1.695   0.0901 .
## exclosure1  -0.29822    0.22657  -1.316   0.1881  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[43]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    161.1    175.0    -72.5    145.1       34 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2145 -0.6554 -0.1027  0.6683  2.0972 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.05459  0.2336  
##  plot             (Intercept) 0.00000  0.0000  
## Number of obs: 42, groups:  quad.unique:plot, 40; plot, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.74580    0.32590   2.288   0.0221 *
## census16su  -0.53419    0.36870  -1.449   0.1474  
## treatmentF   0.01234    0.40387   0.030   0.9756  
## treatmentI   0.52210    0.33381   1.564   0.1178  
## treatmentC   0.05081    0.35508   0.143   0.8862  
## exclosure1  -0.10303    0.22975  -0.448   0.6538  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[44]]
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: no.ind. ~ census + treatment + exclosure + (1 | plot/quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    497.2    527.6   -238.6    477.2      145 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1256 -0.5191 -0.2870  0.2742  2.2798 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  quad.unique:plot (Intercept) 0.027415 0.16558 
##  plot             (Intercept) 0.005977 0.07731 
## Number of obs: 155, groups:  quad.unique:plot, 78; plot, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.52256    0.18194   2.872  0.00408 **
## census16fa  -0.05325    0.17892  -0.298  0.76600   
## census16sp  -0.32291    0.26322  -1.227  0.21991   
## census16su   0.44374    0.15560   2.852  0.00435 **
## treatmentF   0.01956    0.17475   0.112  0.91089   
## treatmentI  -0.21008    0.19612  -1.071  0.28407   
## treatmentC   0.02537    0.18650   0.136  0.89179   
## exclosure1  -0.21845    0.17335  -1.260  0.20760   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Overall, pesticide treatments didn’t affect diversity but abundance. However, when density was tested in individual species, pesticide caused some effects. Snow removal treatment changed both abundance and diversity. I haven’t done the analyses of snow removal on the density of individual species. As to the analysis of pestcide treatment on density of individual species, results were complicated but interesting.